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ABSTRACT 

Stars can produce steady-state winds through radiative driving as long as the mechan- 
ical luminosity of the wind does not exceed the radiative luminosity at its base. This 
upper bound on the mass loss rate is known as the photon-tiring limit. Once above 
this limit, the radiation field is unable to lift all the material out of the gravitational 
potential of the star, such that only part of it can escape and reach infinity. The rest 
stalls and falls back toward the stellar surface, making a steady-state wind impos- 
sible. Photon-tiring is not an issue for line-driven winds since they cannot achieve 
sufficiently high mass loss rates. It can however become important if the star ex- 
ceeds the Eddington limit and continuum interaction becomes the dominant driving 
mechanism. 

This paper investigates the time-dependent behaviour of stellar winds that ex- 
ceed the photon-tiring limit, using 1-D numerical simulations of a porosity moder- 
ated, continuum-driven stellar wind. We find that the regions close to the star show 
a hierarchical pattern of high density shells moving back and forth, unable to escape 
the gravitational potential of the star. At larger distances, the flow eventually be- 
comes uniformly outward, though still quite variable. Typically, these winds have a 
very high density but a terminal flow speed well below the escape speed at the stellar 
surface. Since most of the radiative luminosity of the star is used to drive the stellar 
wind, such stars would appear much dimmer than expected from the super-Eddington 
energy generation at their core. The visible luminosity typically constitutes less then 
half of the total energy flow and can become as low as ten percent or less for those 
stars that exceed the photon-tiring limit by a large margin. 
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1 INTRODUCTION 



The high luminosity of hot, massive stars leads to strong radiatively driven mass loss. In the more 
or less steady winds of O, B, and WR stars, the driving is through scatt ering of radiation by an 
ensemble of spectral lines (e.g., Castor, Abbott, & Klein 1975, hereafter |CAKJ). The tendency of 
these driving lines to become saturated with increasing density acts as an effective wind regulator, 
limiting the wind mass loss rates to ca. 10~ 4 M©/yr in even the most luminous stars. While poten- 
tially quite important for the star's evolution, this level of mass loss is energetically of only minor 
significance, with a total (kinetic and potential) mechanical luminosity of only a few percent of 
the "photon tiring limit" set by the total available stellar luminosity. 

By contrast, the mass loss rates during the giant eruptive phase of Luminous Blue Variable 
(LBV) stars like r] Carinae are inferred to be as high as 0.1 — 1 M^/yr, much higher than can 



be reasonably explained by models based on line-driving (ISmith & Owocki 



2006). In the case of 



r] Carinae, the inferred average mechanical luminosity during its ca. decade-long eruption in the 
1840's is comparable to the estimated average radiative luminosity, ~ 2.5 x 10 7 L e , during this 
time. This luminosity likely well exceeds the Eddington limit, at which the radiative force asso- 
ciated with continuum scattering by free electrons exceeds the inward force of gravity, implying 
then that the extreme mass loss can be driven by continuum, instead of line, opacity. 

Unlike line driving, continuum driving does not suffer saturation at high densities, meaning 
then it can initiate an arbitrarily large mass flux, since all available photon-energy can be used to 
drive the matter. However, if the mechanical luminosity needed to lift this mass flux fully out of 
the gravitational potential exceeds the stellar luminosity, then the "tiring" or reduction in radiative 
energy flux will make it unable to sustain the outflow. This paper presents time-dependent hydro- 
dynamical simulations of the complex combination of stagnated outflow and subsequent infall that 
occurs in continuum-driven models with initial, base flux that exceeds this photon-tiring limit. 

Previous analyses of continuum-driven mass loss have focussed on developing steady-state 
wind outflow models with mass loss rates that approach, but do not exceed the photon-tiring limit. 



To provide the necessary driving regulation and limitation of the mass flux , such mode 



a "porosity" moderation of the opacity associated with a clumped medium (IShav iv 1998 



s invoke 



2001a) 



Such clumps can result from instabilities, such as the photon-bubble instability (ISpiegel 



19761 : 
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Hsu. Arons & Klein 



19971 : 



Spiegel & Tao 



1999; 



Begelman 



2001 



Shaviv 



2001b) in the stellar at- 



mosphere and the upper layers of the stellar envelope. These "photon-bubbles" are a phenomenon 
that occurs in super-Eddington atmospheres and can support large density inhomogeneities. 

This reduces the continuum driving in the dense layers where clumps are optically thick, thus 
allowing the stellar base to be hydrostatically bound even though the luminosity is formally above 
the Eddington limit. But as the clumps become optically thin in the higher-level, lower-density 
layers, they are exposed to the full radiative driving of the assumed super-Eddington luminosity. 
The intermediate layer where the net outward force equals gravity then becomes the sonic point of 
a wind outflow, with the density at this point setting the wind mass loss rate. This identification of 
the sonic radius can be used to predi ct mass loss rates that are consistent with the observations of 



clear-cut super-Eddington objects (IShaviv 



2001a ). 



For a given Eddington parameter T, the mass loss initiated in such models depends on the 
nature of the assumed porosity. Specifically, within the "power-law porosity" formalism developed 



by Owocki, Gayley & Shaviv (2004; hereafter 



OGSJ), it depends on a power- index and the ratio 



of a characteristic "porosity l ength " to the local scale height (see 



OGS 



Analytic solutions derived by 



OGS 



and section 2.4 below). 



show that physically plausible porosity parameters combined 
with moderate super-Eddington parameters, F > 1 can indeed lead to mass loss rates comparable 
to that inferred for the extreme LBV giant eruption of rj Carinae, even approaching the "photon- 
tiring" limit. 



A recent companion paper (van Marie, Owocki & Shaviv 2008; hereafter 



Paper lr ) used numer- 



ic al radi ation hydrodynamics simulations to test these analytic, porosity-moderated wind solutions 



of 



OGS 



for cases that do not exceed this photon-tiring limit. The results showed that the asymp- 
totic steady-states of the numerical simulations are generally in quite close agreement with the 
analytic results. However, with only modest variations in the choice of the key porosity parame- 
ters, these analytic scaling laws also allow base mass fluxes that can exceed the photon-tiring limit, 
with the consequence th at the initial wind solution must necessarilly stagnate at some finite radius 



(van Ma rie et al. 



2( )08aL thereby precluding a net. steady mass outflow. 



The central thrust of the current paper is to explore the nature of the necessarily time-dependent 
flow that ensues from such a porosity model with a base mass flux above this tiring limit. In prin- 
ciple this should involve complex, three-dimensional (3-D) patterns of both outflow and infall. 
But given the computational and conceptual challenges in developing dynamically realistic 3-D 
description, we first explore here much simpler, albeit idealized one-dimensional (1-D) models. A 
key rationale for this approach is that the underlying phenomenon of photon-tiring is not inher- 
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ently multi-dimensional, but rather depends largely on the radial variation of driving vs. gravity. 
The insight gained in studying 1-D simulations of the complex, time-dependent combination of 
outflowing and inflowing layers will form a basis for future, more complete models that include 
multi-dimensional (2-D or even 3-D) flow. 

As detailed in §2 (and also the appendix), a key feature of our simplified method here is to 
maintain a global energy conservation that accounts both for the loss of radiative luminosity in 
regions of outward acceleration, as well as the re-energization of the luminosity from regions of 
infall and shock heating. The results in §3 illustrate the complex pattern of inflow and outflow, 
including also the strong effect on the emergent radiative luminosity. We follow in §4 with a 
discussion of what these results mean for understanding the nature of LBVs, and then conclude in 
§5 with a summary of the work thus far and a brief discussion of directions for future work. 



2 METHOD 

2.1 Mass and momentum conservation equations 

The spherical mass loss simulations presented here are based on evolving the 1-D, time-dependent 
equations for conservation of mass and momentum. For density p and radial flow speed v, mass 
conservation in radius r and time t requires 

Dp _ dp ^ dp p d(vr 2 ) 

Dt dt dr r 2 dr ' 
where D/Dt = d/dt + v d/dr is the total time derivative following the flow. The total change in 

speed depends in turn on the net force-per-unit mass from gas pressure P, gravity, and radiation 

Dv dv dv 1 dP GM 
Dt at or p or r 1 
where, following standard notation, G and M are the gravitation constant and stellar mass, and 

g ra d is the radiative acceleration, described further in §2.4. 

The pressure obeys the ideal gas equation of state, 

P = p— = pa 2 , (3) 
p 

with the latter equality defining the isothermal sound speed a in terms of the temperature T, molec- 
ular weight p, and Boltzmann's constant k. Because a <C t> csc , where v csc = ^2GM/R is the 
escape speed from the stellar surface radius R, gas pressure generally plays little direct role in the 
outward driving of a stellar wind mass loss. But it does provide a basis for matching the outflow 
onto a hydrostatically stratified atmosphere at the wind base. Moreover, as discussed below (§2.2), 
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when an outflow stalls or stagnates, gas pressure plays a key role in dissipating and stopping the 
subsequent gravitational infall. 

As such, our simulations here retain a non-zero gas pressure, but to avoid the necessity of 
modeling in detail the evolution of the (relatively insignificant) gas internal energy, we simply 
assume the gas temperature is a fixed constant. This is set to give a sound speed a = 20 km/s, 
implying a temperature roughly equal to the stellar effective temperature, T ps T eff « 50,000 K. 

2.2 Photon tiring formulation 

A key feature of the present simulations is to account properly for the work done by radiation in 
accelerating material outward against gravity. As detailed in the Appendix, a full treatment should 
in general include the associated time-dependent equations for conservation of radiative momen- 
tum and energy [eqns. (IA4I) and (IA5I) 1. as well as a description of the energy exchange between 
gas and radiation [eqn. (IA3I) 1. However, the further analysis there shows that, by assuming a limit 
of infinite light speed in such a way that both the propagation and diffusion times of radiation 
become small compared to any competing dynamical scale, the storage of energy in the radiation 
energy density E can be effectively bypassed in favor of direct, essentially instantaneous changes 
in the radiative flux F. 

This makes it possible to account for the net work done by the radiation in terms of a simple 
ordinary differential equation (ODE) for the radial change of the radiative luminosity L = 4%r 2 F 
[eqn. <EG5], 

^ = -M g rad - 47rr 2 Q , (4) 
ar 

where the local mass flux, M = Arrpyr 2 , is not necessarily constant (or even positive) in a time- 
dependent flow. The first term on the right side represents the change in luminosity associated with 
the net work done by the radiative acceleration, while the second term represents the luminosity 
change associated with the volume heat exchange rate Q between radiation and the gas energy 
[see eqn. (|A9tl. 

In regions of outflow, the local mass flux is positive, and so the first term represents the basic 
"tiring" effect that leads to a radial decline in the radiative luminosity L. Likewise, in such regions 
of outflow, the need to keep the gas temperature (and internal energy e) constant against the ex- 
pansion cooling requires Q > O [see eqn. (|A9R implying that the second term also contributes to 
a (relatively small) further radial reduction in luminosity. 

In a time-dependent flow, however, the mass flux is not generally constant. Indeed, in regions 
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of flow stagnation and inflow, where v < 0, it can even be locally negative, and therefore lead to a 
radial increase, or re-energization, of the luminosity in the regions above a local infall. Moreover, 
such downflows are typically stopped eventually by strong, localized shock compressions, which 
effectively convert the infall kinetic energy back into radiation, making Q < 0. This represents 
yet another mechanism for re-energization of the luminosity, now occuring at radii above a shock 
compression. 

Using eqn. (IA9I) for the heating rate needed to keep the gas isothermal, the luminosity change 
can be alternatively written as 

dL 9 d(vr 2 ) 

— = -47ipvr 2 g rad - 4vrP^— ^ . (5) 
dr dr 

This equation provides a convenient form to compute the radial change of luminosity in a spher- 
ically symmetric, isothermal flow, applicable even in the fully time-dependent cases considered 
below. At each time- step of the simulation, this ODE can be integrated numerically from the base 
radius R to give the radial run of luminosity, L(r), which is then used to compute the radiative ac- 
celeration g ra d for the next step of the simulation. To avoid the unphysical possibility of a negative 
luminosity, the integration of eqn. © sets a floor luminosity L = 0. 

In practice, maintaining a global energy conservation when there are regions and intervals of 
infall requires one to account also for the energy associated with material that falls back onto the 
stellar surface through the lower boundary radius R. We minimize this here by limiting the inflow 
velocity at the boundary to the sound speed, which then generally leads to shock formation and 
recapture of most of the infall kinetic energy through the instantaneous emission terms described 
above. Nonetheless, there is still a small energy loss through the boundary that has to be accounted 
for to maintain global energy conservation. We do this here by accumulating a storage energy E ac 
that is then assumed to be re-emitted on a parameterized time scale t em , taken here to be several 
(viz. 10) wind flow times. Based on the value of the stored energy at any given time step, the lower 
boundary condition for the luminosity integration of eqn. © then becomes 

L(R) =L* + f^. (6) 



Generally, the associated modification in base luminosity is small, no more than a few percent. 
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2.3 Tiring in steady-state outflows 

For a simple, steady-state flow with constant mass flux M, we can rewrite the radial change of 
luminosity in the forms, 



dL 

dr 



-M 



£rad + P 



d fl 



(GM dv\ 



dr \p 



(7) 



(8) 



where the latter equality assumes an isothermal flow and uses the momentum equation © to 
eliminate the radiative acceleration, with the result that the net tiring no longer depends explicitly 
on the pressure work term. Upon integration from the base radius R, where flow speed is negligibly 
small (v(R) <C a) and the luminosity is set to the interior value, L(R) = L*, the radial dependence 
of luminosity is given by 



L(r) = L* - M 



GM 



1 



R 

r 



(9) 



2 R 

Note in particular that, in this approximation of isothermal flow with constant internal gas energy, 
the pressure work term makes no net contribution to the change in luminosity. 

The form © is in fact the form used in the steady-state, continuum-driven models derived in 



OGSl It also helps justify neglecting the pressure work terms in the time-dependent simulations of 
paper 1, which were used to relax towards the asymptotic steady-states. 
For a wind with terminal speed v^, the emergent luminosity is 



L*-M 



vL GM 
2 R 



(10) 



Setting L c 



the n defines a photon-tiring 



for a steady- state outflow (lOwocki & Gayley 



1997 



imit for the maximum possible mass loss rate 



OGS), 



M tir = 



- — • (11) 

GM/R V ; 

Any flow with base mass flux greater than this limit must necessarily stagnate at a finite radius. 
The simulations presented below examine the complex, time-dependent mixture between outflow 
and infall that results in flows above this tiring limit. 



2.4 Porosity-moderated continuum driving 

The treatment of radiative driving in the time-dependent simuations here is based on th e form alism 



for porosity-moderated, continuum driving developed for the steady-state analysis of 

local luminosity L(r), the radiative acceleration is given by 

K eS L(r) 



OGS 



47J-J-2 c 



For a 



(12) 
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where c is the speed of light, and k cS is an effective opacity that is reduced from the microscopic 
value k by the porous nature of the medium. For a power-law distribution of clump strengths 
characterized by a power index a and a "porosity length" h Q that sets a maximum clump o ptical 



thickness r Q = Kph = pj p Q , the opacity reduction factor depends on the density through Icf. IPGS 
eqn. (57)], 

Hp) = — = — T\ w • ( 13) 

k (1 - a)p/p 

Note that for low densities p <C p a , the clumps all become optically thin, making k — * 1, since the 
material is exposed to the full opacity k. However, at very high densities p ^> p a , self-shadowing of 
optically thick clumps reduces the effective opacity, following an overall power-law dependence 
on density, k ~ (p Q /p) a . The self- shadowing of optically thick clumps is inherently a multi- 
dimensional effect at the microscopic level, but at the microscopic level appears as ID. 

Even in a star for which the continuum opacity k implies an Eddington parameter T = 
kL/A-kGMc above the Eddington limit, i.e. T > 1, the porosity reduction of the radiative force can 
allow high-density regions to remain gravitationally bound, while the outer, lower-density regions 
are driven into a wind outflow. The transition occurs at the sonic point, where r&[p s ] = 1, thus 
defining a base mass loss rate M = An p^aR 2 . For the canonical case a = 1/2 assumed in the 



simulations here, this has the scaling [cf. lOGSl eqn. (77)], 



M = 4(r-1) — , (14) 
r] ac 

where 77 = h a /H, with H = a 2 R 2 /GM being the density scale height at the subsonic base of the 



wind. The ratio of this mass loss rate to the tiring limit given in eqn. (fTTT) has the scaling [cf. 
eqn. (78)], 



OGS 



M T - 1 M/M & 
mtn = -r- = 0.13 (15) 

M tir V a 2o R/Rq 



where a 2 o = a/ 20 km/s. 



2.5 Parameters and numerical implementation 

The simulations here are all based on a single representative stellar model with radius R = 5Oi? , 
and mass M = 50R Q , implying a solar value for the escape speed, v csc = 617 km s -1 . As noted, 
the sound speed is fixed at a = 20 km s _1 , corresponding roughly to a temperature of ca. 50,000 K. 

A summary of the varied model parameters is given in Table 1. Models A-C assume a stellar 
luminosity L* = 1.6 x 1O 7 L , which with an assumed opacity k = 0.4cm 2 /g, implies a strongly 
super-Eddington condition with T = 10. Model D assume a moderate T = 5, while Models E & 
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F assume a marginally super-Eddington condition T = 2. The models are further distinguished 
by the choices of the porosity length parameter r\, as specified in third column of Table 1. For the 
assumed porosity power index a p = 1/2, eqn. (fT5l) shows that the values of 77 and T combine to 
yield a predicted mass-loss tiring ratio m tir , as given in the fourth column of Table 1. As shown 
below, the level by which this base mass-loss tiring exceeds unity is key to the overall nature of 
the flow stagnation and infall of the various models. 

For all cases, we set the re-release timescale for the inner-surface infall energy to be 10 stellar- 
freefall times, t cm = Rj v esc ~ 5.5 x 10 5 s. Since the energy storage and relative luminosity 
contribution is generally small compared to the stellar luminosity, we find the results are generally 
insensitive to the specific release time chosen. 

As in Paper 1 , the simulations here use the ZEUS 3D code developed by 



Stone & Norman 



1992) and 



Clark ( 



1996). But unlike the purely isothermal models previously relaxed to an asymp- 
totic steady-state, the complex time-dependent evolution of photon-tired flows here requires keep- 
ing careful account of the energy exchange, Q, between the gas and radiation. Our implementation 
of ZEUS achieves this by advancing each individual time step assuming adiabatic conditions, with 
an adiabatic index 7 = 5/3, but then resetting the temperature at each grid point back to the as- 
sumed constant value. The difference in internal energy between the constant temperature and the 
temperature that follows from the adiabatic calculations is then added to or subtracted from the 
local radiative luminosity, following eqn. ©. 

In wind simulations, a key general issue is to resolve the acceleration from a nearly hydrostatic 
layer supported by gas pressure, including for example, at the subsonic wind base, which has 
a characteristic scale height H w (2a 2 / v esc 2 )R = 0.002i?. However, in the models here, the 
repeated flow stagnation can often lead to multiple such static layers distributed over many stellar 
radii from the wind base. Moreover, the ultimate escape from the system is only ensured once the 
outward flow speed exceeds the local escape speed, which can sometimes extend to over 100 R. 
To properly resolve this complex flow over a large spatial range, we use n r = 20, 000 spatial zones 
extending over r = 1 — 200 R. To concentrate resolution more in the inner wind, we use a spatial 
grid size that increases by a factor 1 .0005 per zone, i.e. aVj = aVj_i * 1.0005. This leads to a spatial 
resolution that ranges from dr m 5 x 10~ 6 R at r = 1 R to dr nr = 0.1 R at r = 200 R. 

Finally, as in Paper I, the simulations start with an initial condition set by a simple /3-law 
velocity profile, 

r \P 



v(r) = Voo (l - ^ 
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Figure 1. Snapshots of the velocity (a) and log density (b) vs. radius at multiples of 10 5 s after the switching on of photon tiring, at t = 0, for 
Model B. The model was first relaxed to a steady state ignoring tiring. The reduction in the radiative driving flux induces a slow-down towards flow 
stagnation, with shock collisions leading to a build up of high-density shells that eventually fall back onto the star. 



with (3 = 0.5, and a mass flux defined by the porosit y length formalism [eqn. (|14l) l. 

The boundary conditions are the same as used in Paper ll . The density at the inner boundary is 
fixed to a value about 10 times the sonic point density for the analytic porosity model. The velocity 
itself is allowed to float, requiring only that the velocity gradient remain constant across the inner 
boundary, and that the absolute value of the base velocity is limited to the sound speed. Such 



boundary co nditions have been used succesfu lly in simulations of 



in both one (|Owocki. Castor & Rybicki 



ine-driven winds simulations 



1988|) and two dimensions (|Owocki. Cranmer & Blondin 



1994). 



3 RESULTS 

3.1 Initial relaxation and emergence of photon- tiring 

As a first illustration of the effect of photon tiring, figs. [Qplot the radial variation of velocity (a) 
and log-density (b) at fixed time snapshots after turning on the tiring effect in a model which was 
first relaxed to a steady-state without tiring. The model parameters are those of model B, with a 
base mass flux m = 2.34 times the tiring limit. Shortly after the tiring effect is implemented, the 
existing flow begins to decelerate, leading to stagnation and then infall back toward the star. As 
this infalling material meets the remaining upflow below, strong shocks form, compressing the 
material into dense shells, marked by the density spikes in fig. [Tp. Since these dense shells have 
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time |yr] time |yr] 



Figure 2. Left: Colorscale showing mass flux vs. radius and time for Model B. Mass leaves the surface layers only to fall back. Part of the mass 
reaches a higher altitude, where the process repeats itself several times, until the gas can coast away from the star and eventually reach the local 
escape speed. Right: Colorscale showing the associated radiative luminosity as a function of time and radius. The superposed lines depict the 
position of individual mass elements (constant Lagrangian mass). The outward motion of the mass elements increases whenever the luminosity is 
high and slows down, or even reverses if the luminosity is low. But eventually, mass parcels at large radii move systematically outward. 

too much mass to be driven outward by the available radiative flux, they are pulled back by gravity, 
and re-accreted back onto the star. 

Eventually, some of the energy generated can lead to a net outward driving of some of the gas, 
but it too can stagnate before being able to escape, forming more shock layers, and so forth. The 
actual details are of course complex and depend on the specifics of the overloaded initial condition. 

However, rather than dwelling on the specific details which depend on the initial readjustment, 
we shall proceed to focus on the steady-state flow properties, long after the wind relaxes from its 
initial conditions. 

The timescale over which this relaxed state is reached is on the order of several times 10 8 
seconds. This is effectively the crossing time for gas elements with a characteristic velocity of 
approx. 50 km s _1 (the flow speed at the outer boundary) over the length of the grid (200 stellar 
radii). 

3.2 Steady-state behaviour of the shock layers 

To illustrate the complex pattern of infall and outflow in the regions with the overloaded wind, 
fig. [2k presents a colorscale plot of the mass flux for model B (base mass flux 2.34 x M t ir) over 
the radial range r = 1 — 15 R and in a representative time interval of about a third of a year, 
long after relaxation from the initial conditions is reached. The figure shows quite vividly the 
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complex hierarchy of initial base outflows (in red) that are truncated by broad downflows (in blue) 
of stagnated material. This flow pattern is a direct consequence of the photon tiring feedback. 

Indeed, the colorscale in fig. Eb reveals that this flow pattern is tightly coupled with corre- 
sponding time and space variations of the local luminosity. The photon tiring in regions of outward 
acceleration leads to a strong reduction in the luminosity (dark shading); but this is punctuated by 
re-energization of this luminosity (light shading) in the shock collision regions where infall energy 
is converted back into radiation. 

To further demonstrate the overall spatial progression of mass in this complex flow pattern, 
the overlaid light contours in fig. [2b mark the radial and temporal evolution of selected mass 
shells. Because there can be no crossing of such mass shells in these idealized, 1-D, purely radial 
flow simulations, they can be effectively tracked by defining the time and radius variation of a 
Lagrangian mass coordinate, 

"Rcc 



o 



mL(r,t) —An r' 2 p{r' ,t)dr' 

J r 

+ 4tt / r 2 p(r,t')v{r,t')dt', (16) 



with being the outer boundary of the grid (|Owocki. Castor & Rybicki 



1988 



Owocki & Puis 



1999J). The first integral labels each mass shell of the initial (t = 0) condition by the amount of 



mass to the outer boundary; the second integral then tracks the amount of mass that has passed 
through a zone with radius r in the time t since that initial condition. The mass shell tracks plotted 
in fig. [2b therefore represent contours of this Lagrangian mass vs. radius and time. 

Together, the plots in fig. [2] provide an interesting complementary perspective for the flow 
patterns in models with a base mass flux above the tiring limit. Matter is launched from the surface 
of the star, only to fall back almost immediately as the radiation field becomes overloaded and 
can no longer support the outflow against gravity. Once, however, the infalling gas is shocked, it 
releases its kinetic energy as radiation, which can then provide a boost to the overlying material. 
This implies that some of the gas can accelerate above the first shock layer. 

Nevertheless, because the gas accelerating above the first layer is too overloaded, it will again 
stagnate and begin to fall back. The energy it releases can then accelerate gas above the second 
shock, and so forth. However, because progressively higher shocks are further from the star and 
involve a smaller mass flux, their time scale for acceleration and stagnation is longer. 

At a high enough layer, both the total mass flux and the gravitational potential depth are small 
enough such that the remaining luminosity can drive the gas to infinity. Beyond this layer, a freely 
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time [yr] 

Figure 3. Mass flux relative to the tiring mass loss rate (jagged black line), wind velocity relative to the escape speed at the stellar surface 
(dashed line) and radiative luminosity relative to the stellar luminosity (continuous line) at the outer boundary radius 200 R* for model B with 
M = 2.34M t i r , after the simulation has fully relaxed. The mass flux varies slightly below to the tiring limit, while the wind velocity is nearly 
constant at a value well below the surface escape speed. The emergent radiative luminosity is highly variable, but with an average value that is only 
about 10% of the base luminosity, i.e., Too ~ 1. 

coasting wind is accelerated to infinity. In model B, for example, this takes place at roughly 8 
stellar radii. 

3.3 Wind behaviour at larger radii 

To demonstrate the nature of this net outflow at larger radii, fig. [3] plots the time variation of 
three key quantities at the outer boundary of the computational grid (200 R*) for the relaxed state 
of model B. These include the radiative luminosity relative to the stellar luminosity, the wind 
velocity relative to the escape speed, and the mass flux relative to the photon-tiring mass loss rate. 
The driving proves to be quite efficient, bringing the mass flux close to (~ 85%) the tiring limit. 
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The wind velocity is nearly constant at ca. 50km s -1 , about a factor 2.5 times the sound speed, 
and much less than the surface escape speed, v csc (R) « 620 km/s; but it is just slightly above the 
local escape speed, v esc (r = 200 R) ~ 44 km/s. The mass flux is also quite constant, with the 
small fluctuations that mostly reflect residual variations in the flow density from the discrete shells 
formed near the wind base. 

However, note that, at least in this idealized ID model, the radiative luminosity exhibits quite 
strong variations. The average level is also subtantially reduced, to about ~ 10% of the base 
luminosity. But there are also brief spikes above 50%, arising from the sudden radiative emissions 
in the shock infall regions. More realistically, such spikes would likely to be smoothed out by both 
radiative storage and 3-D averaging effects ignored in the present 1-D, sudden radiative release 
simulations. 

Finally, fig. @] compares the outer boundary state of models C, B, and A with increasing mass 
flux tiring factors, 1.17, 2.34, and 4.68. To show finer details of the variations, the time resolution 
in fig. [4] is much higher than in fig. [3] Qualitatively, the behaviour of the wind is similar for all 
three simulations: the mass flux is a large fraction (> 65% of the tiring limit in all three cases). 
The wind moves slowly, and and the radiative luminosity that escapes at the outer boundary is 
highly variable but greatly reduced from the base value. 

The results also reveal notable trends in the quantitative properties. In particular, for models 
that exceed the tiring limit by a greater margin, the net driving of mass loss becomes somewhat 
less efficient. More photons escape while the mass flux decreases and becomes less variable. This 
is understandable as the motion of matter back and forth close to the star recycles energy from 
the radiation field, through kinetic and potential energy of the wind, and back into radiation. If 
the wind exceeds the tiring limit by a wider margin, this process become more extensive, which 
reduces the efficiency of the driving and allows more photons to escape. However, the outer wind 
velocity is almost the same in all three models. 

If we compare models with a high T and high tiring number, to models with low T and low 
tiring number, then we see a clear trend. With a smaller tiring number, there are fewer shock layers 
up to where free outflow is reached. As a consequence, the light curve is much more regular, as 
can be expected since it is composed of fewer shock "oscillations". This can be seen in figs. [6] and 
|5l A lower luminosity (lower T) means that a smaller amount of mass is present in the grid at any 
given time. Therefore, the optical depth is smaller, so the photons have a better chance to escape. 
As a result, lower luminosities at the base generally lead to lower efficiency in the driving. 

In Table 1, we show the relative contributions of potential energy luminosity (E g = l/2M(v esc (R) 2 — 
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time [yr] 

Figure 4. The same variables as fig.[3]for models C, B, and A, having (from left to right) M=1.17, 2.34 and 4.68M t i r respectively. As the mass 
flux at the surface exceeds the tiring limit by a larger margin, the tiring mechanism becomes somewhat less efficient: the mass flux at the outer 
boundary (jagged black line) decreases and radiative luminosity (spiky, continuous line) increases. Also, the variations in mass flux smooth out and 
decrease in size. The wind speed (dashed line) remains nearly the same for all three simulations. 



Table 1. Summary of model parameters and results. The other model parameters are M = 50Mq, R = 50Rq, k = 0.4cm 2 gr~ 1 giving 
L E( jd = 1.6 10 b Lq, a = 20 km s~ 1 and a p = 1/2. For the T = 10 models, M t i r = 0.52 Mq /yr. n max is defined as the maximum number 
of radial shock layers, which quantifies the complexity of the radial structure, while rtop is the largest radius to which shocks reach. L/E and 
Eg I E are the contributions of radiative luminosity and gravitational potential energy respectively to the energy outflow at the outer boundary. 
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0.25 


4.68 


0.66 
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2.34 
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0.7 
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0.08 


8 


0.076 


0.916 
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10 


1.0 


1.17 


0.85 


0.5 
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0.08 
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0.059 


0.935 
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0.25 


2.08 


0.68 


1.5 


5 


0.08 
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0.307 


0.689 
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0.06 


2.17 


0.38 


4.4 


4 


0.09 


2.8 


0.466 


0.527 


F 


2 


0.12 


1.08 


0.47 


4.3 


3 


0.09 


3.5 


0.450 


0.544 



v esc (r) 2 ) and radiative (L) luminosity to the total energy outflow (E) at the outer boundary. These 
two make up the bulk of the energy flow, because kinetic energy flow (l/2Mt;^) and the internal 
energy advection (MkT ' /run) together contribute less than 0.5%. Maintaining a fully accurate en- 
ergy conservation in these simulations can be difficult due to the numerical problems that occur 
when resolving a large number of radiative shocks. We achieved an accuracy of more than 90% in 
all simulations with the exception of simulation A, which has the largest number of shocks, and 
for which the energy conservation accuracy was about ~ 75%. Therefore, although the results of 
the simulation seem to follow the trends observed in this section, the quantitative values should be 
viewed with some caution. 



16 A. J. van Marie et al. 



8E+06 - 



6E+06 - 



© 



4E+06 - 



2E+06 - 



r = 2, T] = 0.12 - 

r = 2, T] = 0.06 




time [yr] 



Figure 5. The light curve of two slightly photon-tired cases, with r; 
light-curves are much more regular, and not chaotic. 



0.06 and r\ = 0.12. Both have T = 2. Unlike the higher T models, the 



4 DISCUSSION 

It should be emphasized here that the photon-tiring effects that are the focus of this paper are only 
relevant for the most extreme forms of radiatively driven mass loss. For even the strongest, steady- 
state, line-driven winds from WR stars, the mass loss rates can range up to a few times 10~ 5 M /yr 
and the flow speeds up to Voo ~ 3000 km/s. This then implies a mechanical wind liminosity that is 
less than a few percent of a typical WR luminosity of ca. 10 6 L Q . 

By contrast, the Homunculus nebulae that resulted from the 1840-1850 giant erupt ion of 



77 Car inae is estimated to harbor at least 10 M with expansion speeds on order 600 km/s (ISmith 



20021) . implying a characteristic average mechanical luminosity that is quite comparable to the 
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Figure 6. The radial dependence of the mass flux as a function of time (i.e., similar to fig[T])), for model F, with T = 2 and M/M t i r = 1.08. 
Because of the more moderate conditions, the shock layer structure is much simpler-it includes fewer layers, extends to a smaller radius, and 
significantly less chaotic (as can be seen in fig. [5) 



estimated radiative luminosity 2.5 x 10 7 L & during the epoch. This seems a strong indication that 
photon tiring may have played a key role in the driving of this extensive mass loss. 

But there are some key differences between these properties of the i] Carinae eruption and 
the 1-D photon-tired wind models computed here. For one, the observed terminal flow speed of 
ca. 600 km/s is much higher than the ca. 50 km/s found for the tiring-limited models computed 
here. Moreover, in contrast to the rough parity between mechanical and radiative luminosity for 
T] Carinae's giant eruption, these tiring models give mechanical luminosities that are several times 
higher than the emergent radiation, which itself is also quite variable. 

However, as already noted, it seems quite likely that these general properties of 1-D photon- 
tiring could be quite different in full 3-D case. In particular, the likely lateral variations between 
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regions of infall and outflow could allow a stronger, more regular escape of radiative flux. Indeed it 
might even resemble the extensive, complex 3-D clump structure envisioned in the parameterized 
porosity formalism. The net effect could be a self-regulation of porosity to keep the mass flux 
below the tiring limit, and thus allow the kind of quasi-steady, higher- speed w ind outlfow obtained 



in both the analtyic models of 



Paper 1 



' OGS and the numerical simulations of ] 
The up-flow and down-flow nature of the mass flux close to the star is also reminiscent of 
convective motion, albeit driven by radiative luminosity rather than buoyancy. In the stellar inte- 
rior, convection can effectively transport energy at such a rate that the star can locally exceed the 
Eddington limit (IJoss^ Salpeter & Ostrikenll973|) . But in surface layers it becomes inefficient due 
to the lower density (lOGSl) . Similarly, the shock layers allow for mechanical energy advection 
to keep the photon-tired layer with a reduced luminosity. However, unlike convection where the 
flows are sub-sonic, the photon tired layer is composed of highly super-sonic motion. 

Even in the 1-D framework, the present simulations assume an instantaneous nature for both 
the propagation and diffusion of radiation, and neglect optically thick backwarming effects in 
favor of a simple, constant temperature. Including an optically thick radiative diffusion with a 
finite speed would allow energy to be stored in the radiation field, with a likely smoothing effect 
on both luminosity variation and the resulting motions of mass shells. In 3-D or even 2-D such 
shells should break up into clumps. This could allow more photons to escape between the dense 
clumps, reducing the mass flux and increasing the radiative luminosity at the outer boundary. This 
higher luminosity at large radii could also accelerate the flow to higher terminal speeds. 



5 SUMMARY & FUTURE WORK 

If the base mass flux from a star exceeds the photon-tiring limit, its behaviour is completely dif- 
ferent from that of ordinary radiation driven winds. Rather than producing a steady wind with a 
velocity exceeding the escape speed at the stellar surface, the tiring results in a complex combi- 
nation of infall and outflow in the inner wind, with a terminal velocity in the outer wind that is 
much lower than the surface escape speed. The net mass flux is somewhat below the photon-tiring 
mass-loss limit. Both quantities seem to depend only weakly on the wind parameters at the stel- 
lar surface as set by the porosity length formalism. (N.B. This may change for multi-dimensional 
simulations.) Less than half of the stellar luminosity escapes in the form of radiation, even for 
stars that are only a small fraction above the photon tiring limit. For stars that are further above the 
photon-tiring limit, this fraction drops rapidly to less than 10%, so stars like this have an effective 



Stellar winds above the photon-tiring limit 19 

radiative luminosity that is less than the Eddington limit and could appear to be quite dim to an 
external observer. 



We note in passing that the flow stagnation and a complex combina tion of infal 



has also proposed for the circumstellar environment of some AGB stars (ISoker 



2008). 



and outlfow 



In future work we plan to introduce time-dependent, flux-limited, radiative diffusion in our 
simulations, first in a 1-D model, but eventually with aims to extend this to a multi-dimensional 
grid. 

We also intend to investigate the effect of stellar rotation on the geometry of the continuum- 
driven wind, and how this might give rise to the bipolar shape seen the homonculus nebula of 
T) Carinae. 
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APPENDIX A: REDUCED CONSERVATION EQUATIONS FOR RADIATION AND 

GAS 



Startin g from the genera l mixe d-fram e equations of radiatio n hydrodynamics as given by 
Shaviv ( 2001al) and 



Mihalas & Mihalas 



(1984) [see also 



Spi egel & Tad (|1999|) 1. we first write the basic equations for 



conservation of mass, momentum, and energy for a gas of density p, pressure P, specific internal 
energy e = 3P/2p, and velocity u, 



De 
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P- 
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pV ■ u 

VP 
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srad 
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Dt \p 
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(Al) 
(A2) 

(A3) 



where D/Dt = d/dt + u-V represents the advective time derivative. Here g ra d and g grav are the 
accelerations associated with radiative driving and gravity respectively, while Q is the net rate per 
unit volume for either radiative heating (Q > 0) or cooling (Q < 0). 

Likewise, for radiation with energy density E and energy flux F, the associated conservation 
equations for radiative momentum and energy can be written, 

PkF 



1 dF 1 ^ 

? 1ft + 3 VS 



-Pgrad 



BE 



(A4) 
(A5) 



-Q - pu ■ g rad , 

where eqn. (IA4I) assumes the Eddington approximation in setting the radiation pressure to E/3, 
with the extra equality giving the radiative acceleration in terms of the flux F and a grey opacity 
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k, ignoring radiative drag terms of the form uE/c for the highly subrelativistic (u <C c) flows 
considered here. The right side of eqn. (IA5I) accounts for the radiative energy expended in heating 
and accelerating the gas, the latter representing the key "photon tiring" effect that is the focus of 
this paper. 

In practice, a further approximation here is to ignore the partial time derivatives in both the ra- 
diative momentum and energy equations (IA4|) and (|A5I) . For the radiative flux term, this amounts 
to ignoring light travel signals across the system, effectively assuming the instantaneous propaga- 
tion of the radiative momentum in a limit of infinite light speed, c — > oo. For the radiative energy 
term, it amounts to also assuming an instantaneous diffusion of radiative energy, even compared to 
hydrodynamical variations. For a layer of optical depth r, the ratio of diffusion to dynamical time 
scales scales roughly as ra/c, and since the sound speed is a ~ 20 km/s, this assumption should 
hold for r <C 10 4 . Even for the very dense outflows considered here, this condition is generally 
well satisified. 

In the standard case of a static, planar, grey atmosphere in radiative equilibrium, the right side 
of eqn. (IA5I) vanishes, implying a constant radiative flux and allowing eqn. (IA4I) to be trivially 
integrated to yield a simple optical depth scaling between radiative energy and flux, E ~ tF. The 
exchange of energy with the gas then leads to an equilibriation of the radiative and gas temperature, 
resulting in an effective "backwarming" of the gas at large optical depths. 

An analogous backwarming can also be expected in a very dense, optically thick wind outlow, 
but the overall dynamical effect of this is likely to be quite limited, since the associated gas internal 
energy is still well below the gravitational binding energy. The associated energy ratio can be cast 
in terms of the squared ratio between the sound speed and escape speed, which in a wind with 
optical depth r should scale as 



where a* is the sound speed associated with the star's effective temperature T eS . Since typically 
al/ v esc 2 ~ 1CT 3 , this ratio only becomes of order unity for optical depths r ~ 10 12 , characteristic 
of the very center of an entire star. It is always small for any stellar mass outflow. 

Since this implies that gas internal energy is dynamically unimportant, we ignore such optically 



gas at a nearly constant temperature, chosen here to give a fixed sound speed a = 20 km/s that 
is characteristic of gas near the stellar effective temperature T e g . Since this implies a constant gas 




(A6) 



thick backwarming effects, and instead assume that the radiative heating term Q acts to keep the 
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internal energy e = 3a 2 /2, the energy equations for the gas (IA3I) and radiation (IA5I) reduce to 

- ^7 = -Q = - PV • u (A7) 
V • F = -Q-pu- g rad , (A8) 

where the latter equality in eqn. (IA7I) follows from the mass conservation equation (lAlj) . For the 
present simple case of a 1-D spherical flow with radial speed v and local luminosity L = Anr 2 F, 
these reduce to 

PD, = _PJg) 

^ = -4vrr 2 g - M(? rad , (A10) 
ar 

where M = A-Kpvr 2 is the local mass flux, which is not necessarily constant (or even positive) in 
a time-dependent wind. 

Note that, in this simplified formulation, the radiative momentum equation is superfluous, and 
the radiative internal energy irrelevant. But there is still a clear global conservation of energy 
between the radiation and both the work and the heating it imparts to the gas. As such, this does 
provide a convenient 1-D formalism for examining the nature of radiatively driven mass flows near 
the photon-tiring limit. 



